/*
 * Cluster.cpp
 *
 *  Created on: Jul 28, 2011
 *      Author: nzhao
 */

#include "core/Cluster.h"


void coutHami(double * hami, int dim, int clst_idx);
void cout_C_Hami(CLD ** C_HAMI, int blk, int * dim, int idx, int sign);
void cout_complex_mat(gsl_matrix_complex * mat, int n1, int n2);
double Lorentzian(double omega, double omega0, double gamma);

void calculation_cluster(cCLUSTER cluster, cSPIN * spin_list, int order_i, int idx_j, int count, double * clst_res,  double * irr_clst_res)
{
	int nspin, time_idx, hamiform_idx, idx, k;
	int * p_clst_idx;
	double * timeseq, subval;

	nspin 		= order_i + 1;
	p_clst_idx 	= CLST_INDEX[order_i] + idx_j * nspin;
	hamiform_idx = CLUSTER_2_HamiFORM[order_i][idx_j];

	switch (TASK_MODE)
	{
		/////////////////////////////////////////////////////////////////////////////////////
		//Ensemble CCE
		case 0:
			cluster.initialize(spin_list, nspin, p_clst_idx, PULSE_NUM, count, hm[ hamiform_idx ]);
			cluster.calc_Hamiltonian(  1 );
			cluster.calc_Hamiltonian(  2 );
			cluster.pre_Signal();
			for (time_idx = 0; time_idx <= NTIME; time_idx ++)
			{
				timeseq = TIME_SEQ_INTERVAL + time_idx * INTERVAL_NUM;
				clst_res[time_idx] = cluster.calc_Signal( timeseq );
			}

			for (time_idx = 0; time_idx <= NTIME; time_idx ++)
			{
				subval = 1.;
				for(k = 0; k < sblst[count].sublen; k++)
				{
					idx = sblst[count].sublist[k];// - 1;
					assert( idx < count);
					subval *= RES_BUFFER_IRREDUCE[idx * (NTIME+1) + time_idx];
				}

				irr_clst_res[time_idx] = clst_res[time_idx] / subval;
				if(order_i < CALCULATED_ORDER - 1 || CACHE_Q)
					RES_BUFFER_IRREDUCE[count * (NTIME+1) + time_idx] = irr_clst_res[time_idx];
			}
			cluster.clear();
			break;

		/////////////////////////////////////////////////////////////////////////////////////
		//Singel-Sampel CCE
		case 1:

			cluster.initialize_SGL(spin_list, nspin, p_clst_idx, PULSE_NUM, count, hm[ hamiform_idx ], STATE_J_LIST[count]);
			cluster.calc_Hamiltonian_SGL(  1 );
			cluster.calc_Hamiltonian_SGL(  2 );
			cluster.pre_Signal_SGL();
			for (time_idx = 0; time_idx <= NTIME; time_idx ++)
			{
				timeseq = TIME_SEQ_INTERVAL + time_idx * INTERVAL_NUM;
				clst_res[time_idx] = cluster.calc_Signal_SGL( timeseq );
			}

			for (time_idx = 0; time_idx <= NTIME; time_idx ++)
			{
				subval = 1.;
				for(k = 0; k < sblst[count].sublen; k++)
				{
					idx = sblst[count].sublist[k];// - 1;
					assert( idx < count);
					subval *= RES_BUFFER_IRREDUCE[idx * (NTIME+1) + time_idx];
				}

				irr_clst_res[time_idx] = clst_res[time_idx] / subval;
				if(order_i < CALCULATED_ORDER - 1 || CACHE_Q)
					RES_BUFFER_IRREDUCE[count * (NTIME+1) + time_idx] = irr_clst_res[time_idx];
			}

			cluster.clear_SGL();
			break;

		/////////////////////////////////////////////////////////////////////////////////////
		// correlation function
		case 2:
			cluster.initialize_Correlation(spin_list, nspin, p_clst_idx, PULSE_NUM, count, hm[ hamiform_idx ]);
			cluster.calc_Hamiltonian_Correlation( 1 );
			cluster.calc_Hamiltonian_Correlation( 2 );
			cluster.pre_Signal_Correlation();
			for (time_idx = 0; time_idx <= NTIME; time_idx ++)
			{
				clst_res[time_idx] = cluster.CorrelationFunc(MIN_TIME + time_idx * TIME_STEP);
			}

			for (time_idx = 0; time_idx <= NTIME; time_idx ++)
			{
				subval = 0.;
				for(k = 0; k < sblst[count].sublen; k++)
				{
					idx = sblst[count].sublist[k];// - 1;
					assert( idx < count);
					subval += RES_BUFFER_IRREDUCE[idx * (NTIME+1) + time_idx];
				}

				irr_clst_res[time_idx] = clst_res[time_idx] - subval;
				if(order_i < CALCULATED_ORDER - 1 || CACHE_Q)
					RES_BUFFER_IRREDUCE[count * (NTIME+1) + time_idx] = irr_clst_res[time_idx];
			}
			cluster.clear_Correlation();
			break;
		/////////////////////////////////////////////////////////////////////////////////////
		// Gaussian approximation
		case 3:
			cluster.initialize_Correlation(spin_list, nspin, p_clst_idx, PULSE_NUM, count, hm[ hamiform_idx ]);
			cluster.calc_Hamiltonian_Correlation( 1 );
			cluster.calc_Hamiltonian_Correlation( 2 );
			cluster.pre_Signal_Correlation();
			for (time_idx = 0; time_idx <= NTIME; time_idx ++)
			{
				clst_res[time_idx] = cluster.GaussianEcho(MIN_TIME + time_idx * TIME_STEP);
			}

			for (time_idx = 0; time_idx <= NTIME; time_idx ++)
			{
				subval = 1.;
				for(k = 0; k < sblst[count].sublen; k++)
				{
					idx = sblst[count].sublist[k];// - 1;
					assert( idx < count);
					subval *= RES_BUFFER_IRREDUCE[idx * (NTIME+1) + time_idx];
				}

				irr_clst_res[time_idx] = clst_res[time_idx] / subval;
				if(order_i < CALCULATED_ORDER - 1 || CACHE_Q)
					RES_BUFFER_IRREDUCE[count * (NTIME+1) + time_idx] = irr_clst_res[time_idx];
			}
			cluster.clear_Correlation();
			break;

		/////////////////////////////////////////////////////////////////////////////////////
		// Noise spectrum
		case 4:
			cluster.initialize_Correlation(spin_list, nspin, p_clst_idx, PULSE_NUM, count, hm[ hamiform_idx ]);
			cluster.calc_Hamiltonian_Correlation( 1 );
			cluster.calc_Hamiltonian_Correlation( 2 );
			cluster.pre_Signal_Correlation();
			for (time_idx = 0; time_idx <= NTIME; time_idx ++)
			{
				clst_res[time_idx] = cluster.NoiseSpectrum(MIN_TIME + time_idx * TIME_STEP);
			}

			for (time_idx = 0; time_idx <= NTIME; time_idx ++)
			{
				subval = 0.;
				for(k = 0; k < sblst[count].sublen; k++)
				{
					idx = sblst[count].sublist[k];// - 1;
					assert( idx < count);
					subval += RES_BUFFER_IRREDUCE[idx * (NTIME+1) + time_idx];
				}

				irr_clst_res[time_idx] = clst_res[time_idx] - subval;
				if(order_i < CALCULATED_ORDER - 1 || CACHE_Q)
					RES_BUFFER_IRREDUCE[count * (NTIME+1) + time_idx] = irr_clst_res[time_idx];
			}
			cluster.clear_Correlation();
			break;
		/////////////////////////////////////////////////////////////////////////////////////
		// error
		default:
			cout << "ERROR: Wrong TASK_MODE" << endl;
			assert(0);
			break;
	}
}

void cCLUSTER::initialize(cSPIN * spins, int num_spin, int * spinIdx, int npulse, int total_idx, clusterHamiltonian hamiform)
{
	int i, j, k, idx_i, idx_j, idx_ij;
	assert(num_spin <= MAX_SPIN_IN_CLST);

	spinlist		= spins;
	nspin			= num_spin;
	spin_index		= spinIdx;
	pulse_num		= npulse;
	clst_idx		= total_idx;

	sINT_TENSOR1 dip;

	//	dim = 1;
	dim = hamiform.FullDim();
	for (i = 0; i < nspin; i++)
	{
//		dim *= 2;
		idx_i = spin_index[i];

		hf1x[i] = spinlist[idx_i].hf(1).x;
		hf1y[i] = spinlist[idx_i].hf(1).y;
		hf1z[i] = spinlist[idx_i].hf(1).z;

		hf2x[i] = spinlist[idx_i].hf(2).x;
		hf2y[i] = spinlist[idx_i].hf(2).y;
		hf2z[i] = spinlist[idx_i].hf(2).z;

		for (j = 0; j < nspin; j++)
		{
			idx_j = spin_index[j];
			idx_ij = idx_i * NBATH + idx_j;

			dip = calc_dipole_int(idx_i, idx_j);

			if(IsMediated)
			{
			xx1 [i * nspin + j] = dip.xx + MEDIATED_INT1[idx_ij + 0];
			yy1 [i * nspin + j] = dip.yy + MEDIATED_INT1[idx_ij + 1];
			zz1 [i * nspin + j] = dip.zz + MEDIATED_INT1[idx_ij + 2];

			xy1 [i * nspin + j] = dip.xy + MEDIATED_INT1[idx_ij + 3];
			yz1 [i * nspin + j] = dip.yz + MEDIATED_INT1[idx_ij + 4];
			zx1 [i * nspin + j] = dip.zx + MEDIATED_INT1[idx_ij + 5];

			yx1 [i * nspin + j] = dip.yx + MEDIATED_INT1[idx_ij + 6];
			zy1 [i * nspin + j] = dip.zy + MEDIATED_INT1[idx_ij + 7];
			xz1 [i * nspin + j] = dip.xz + MEDIATED_INT1[idx_ij + 8];

			// xx1==xx2...
			xx2 [i * nspin + j] = dip.xx + MEDIATED_INT2[idx_ij + 0];
			yy2 [i * nspin + j] = dip.yy + MEDIATED_INT2[idx_ij + 1];
			zz2 [i * nspin + j] = dip.zz + MEDIATED_INT2[idx_ij + 2];

			xy2 [i * nspin + j] = dip.xy + MEDIATED_INT2[idx_ij + 3];
			yz2 [i * nspin + j] = dip.yz + MEDIATED_INT2[idx_ij + 4];
			zx2 [i * nspin + j] = dip.zx + MEDIATED_INT2[idx_ij + 5];

			yx2 [i * nspin + j] = dip.yx + MEDIATED_INT2[idx_ij + 6];
			zy2 [i * nspin + j] = dip.zy + MEDIATED_INT2[idx_ij + 7];
			xz2 [i * nspin + j] = dip.xz + MEDIATED_INT2[idx_ij + 8];
			}
			else
			{
				xx1 [i * nspin + j] = dip.xx;
				yy1 [i * nspin + j] = dip.yy;
				zz1 [i * nspin + j] = dip.zz;

				xy1 [i * nspin + j] = dip.xy;
				yz1 [i * nspin + j] = dip.yz;
				zx1 [i * nspin + j] = dip.zx;

				yx1 [i * nspin + j] = dip.yx;
				zy1 [i * nspin + j] = dip.zy;
				xz1 [i * nspin + j] = dip.xz;

				// xx1==xx2...
				xx2 [i * nspin + j] = dip.xx;
				yy2 [i * nspin + j] = dip.yy;
				zz2 [i * nspin + j] = dip.zz;

				xy2 [i * nspin + j] = dip.xy;
				yz2 [i * nspin + j] = dip.yz;
				zx2 [i * nspin + j] = dip.zx;

				yx2 [i * nspin + j] = dip.yx;
				zy2 [i * nspin + j] = dip.zy;
				xz2 [i * nspin + j] = dip.xz;
			}
		}
	}
	ZEROS[0] = 0.0;

	set_localfield();

	hm          = hamiform;
	blk_num     = hm.BlockNumer();
	blk_dim     = hm.HamiDim();

	C_HAMI_0 = new CLD * [blk_num];     eigen_val0 = new double * [blk_num];    eigen_vec0 = new CLD * [blk_num];
	C_HAMI_1 = new CLD * [blk_num];     eigen_val1 = new double * [blk_num];    eigen_vec1 = new CLD * [blk_num];
	C_HAMI_2 = new CLD * [blk_num];     eigen_val2 = new double * [blk_num];    eigen_vec2 = new CLD * [blk_num];
	W_matrix = new CLD * [blk_num];

	left_mat0  = new CLD * [blk_num];
	left_mat   = new CLD * [blk_num];
	right_mat0 = new CLD * [blk_num];
	right_mat  = new CLD * [blk_num];

	int n_dim1, n_dim2, n_dim = 0;
	for(i = 0; i < blk_num; i++)
	{
		n_dim += blk_dim[i];
		n_dim1 = blk_dim[i];
		n_dim2 = n_dim1 * n_dim1;
		C_HAMI_0[i] = new CLD [ n_dim2 ];   eigen_val0[i] = new double [n_dim1];    eigen_vec0[i] = new CLD [n_dim2];
		C_HAMI_1[i] = new CLD [ n_dim2 ];   eigen_val1[i] = new double [n_dim1];    eigen_vec1[i] = new CLD [n_dim2];
		C_HAMI_2[i] = new CLD [ n_dim2 ];   eigen_val2[i] = new double [n_dim1];    eigen_vec2[i] = new CLD [n_dim2];
		W_matrix[i] = new CLD [ n_dim2 ];

		left_mat0[i]  = new CLD [n_dim2];
		right_mat0[i] = new CLD [n_dim2];
		left_mat[i]   = new CLD [n_dim2];
		right_mat[i]  = new CLD [n_dim2];
	}
	assert(n_dim == dim);

	i_phase_fact = new CLD [ dim * 2 * (pulse_num + 1)];
	exp_eigen_val = new CLD [dim * 2 * (pulse_num + 1)];
}

void cCLUSTER::clear()
{
    int i;
    for(i = 0; i < blk_num; i++)
    {
        delete [] C_HAMI_0[i];  delete [] eigen_val0[i];    delete [] eigen_vec0[i];
        delete [] C_HAMI_1[i];  delete [] eigen_val1[i];    delete [] eigen_vec1[i];
        delete [] C_HAMI_2[i];  delete [] eigen_val2[i];    delete [] eigen_vec2[i];
        delete [] W_matrix[i];
        delete [] left_mat0[i];		delete [] left_mat[i];
        delete [] right_mat0[i];	delete [] right_mat[i];
    }
    delete [] C_HAMI_0; C_HAMI_0 = NULL;
    delete [] C_HAMI_1; C_HAMI_1 = NULL;
    delete [] C_HAMI_2; C_HAMI_2 = NULL;

    delete [] eigen_val0; eigen_val0 = NULL;
    delete [] eigen_val1; eigen_val1 = NULL;
    delete [] eigen_val2; eigen_val2 = NULL;

    delete [] eigen_vec0; eigen_vec0 = NULL;
    delete [] eigen_vec1; eigen_vec1 = NULL;
    delete [] eigen_vec2; eigen_vec2 = NULL;

    delete [] W_matrix; W_matrix =NULL;
    delete [] left_mat0; left_mat0 = NULL;
    delete [] right_mat0; right_mat0 = NULL;
    delete [] left_mat; left_mat = NULL;
    delete [] right_mat; right_mat = NULL;

    delete [] i_phase_fact; i_phase_fact = NULL;
    delete [] exp_eigen_val; exp_eigen_val = NULL;
}

void cCLUSTER::calc_Hamiltonian(int sign)
{
	int i, j;

        CLD     ** C_HAMI;
        double  ** eigenval;
	CLD 	** eigenvec;

	double * xx, * yy, * zz, * xy, * yz, * zx, * yx, * zy, * xz;

	if(sign == 1)
	{
		for (i = 0; i < nspin; i++)
		{
			Bx[i] = locX[i] + hf1x[i];
			By[i] = locY[i] + hf1y[i];
			Bz[i] = locZ[i] + hf1z[i];

		}
		xx = xx1; yy = yy1; zz = zz1;
		xy = xy1; yz = yz1; zx = zx1;
		yx = yx1; zy = zy1; xz = xz1;

		C_HAMI = C_HAMI_1;
		eigenval = eigen_val1;
		eigenvec = eigen_vec1;
	}
	else if(sign == 2)
	{
		for (i = 0; i < nspin; i++)
		{
			Bx[i] = locX[i] + hf2x[i];
			By[i] = locY[i] + hf2y[i];
			Bz[i] = locZ[i] + hf2z[i];
		}
		xx = xx2; yy = yy2; zz = zz2;
		xy = xy2; yz = yz2; zx = zx2;
		yx = yx2; zy = zy2; xz = xz2;

		C_HAMI = C_HAMI_2;
		eigenval = eigen_val2;
		eigenvec = eigen_vec2;

	}
	else // sign==3
	{
		for (i = 0; i < nspin; i++)
		{
			Bx[i] = locX[i] + 0.5 * (hf1x[i] + hf2x[i]);
			By[i] = locY[i] + 0.5 * (hf1y[i] + hf2y[i]);
			Bz[i] = locZ[i] + 0.5 * (hf1z[i] + hf2z[i]);
		}
		xx = xx2; yy = yy2; zz = zz2;
		xy = xy2; yz = yz2; zx = zx2;
		yx = yx2; zy = zy2; xz = xz2;

		C_HAMI = C_HAMI_0;
		eigenval = eigen_val0;
		eigenvec = eigen_vec0;
	}

        INTERACTIONS[0] = ZEROS;
        INTERACTIONS[1] = Bx; INTERACTIONS[2] = By; INTERACTIONS[3] = Bz;
        INTERACTIONS[4] = xx; INTERACTIONS[5] = yy; INTERACTIONS[6] = zz;
        INTERACTIONS[7] = xy; INTERACTIONS[8] = yz; INTERACTIONS[9] = zx;
    	INTERACTIONS[10] = yx; INTERACTIONS[11] = zy; INTERACTIONS[12] = xz;

        double re, im;
        for(i = 0; i < blk_num; i++)
        {
            for(j = 0; j < blk_dim[i]*blk_dim[i]; j++)
            {
                re = hm.ElementVal_R(INTERACTIONS, i, j);
                im = hm.ElementVal_I(INTERACTIONS, i, j);
                C_HAMI[i][j] = CLD(re, im);
            }
        }
//        if(dim == 4)
//        	cout_C_Hami(C_HAMI, blk_num, blk_dim, clst_idx, sign);

        for(i = 0; i < blk_num; i++)
            MKL_hermite_mat_diag(blk_dim[i], C_HAMI[i], eigenval[i],eigenvec[i]);
}

void cCLUSTER::pre_Signal()
{
    int i;
    char * transa, * transb;

    transa = "ConjugateTranspose";
    transb = "NoTranspose";

    for(i = 0; i < blk_num; i++)
    {
        MKL_mat_mult(blk_dim[i], transa, transb, COMPLEX_ONE, eigen_vec2[i], eigen_vec1[i], COMPLEX_ZERO, W_matrix[i]);
        MKL_matrix_copy('C', blk_dim[i], COMPLEX_ONE, eigen_vec1[i], right_mat0[i]);
        MKL_matrix_copy('N', blk_dim[i], COMPLEX_ONE, eigen_vec2[i], left_mat0[i]);
    }
}

double cCLUSTER::calc_Signal(double * timeseq)
{
	int i, j;
	double	res;

	if(PULSE_SYMMETRY_MODE[0] > 0)
	{
		if(pulse_num % 2 == 0)
		{
			res = PDD_algorithm_even(timeseq);
		}
		else
		{
			res = PDD_algorithm_odd(timeseq);
		}
	}
	else// PULSE_SYMMETRY_MODE[0] == 0
	{
		res = Generic_algorithm(timeseq);
	}

	return res;
}

double cCLUSTER::PDD_algorithm_even(double * timeseq)
{
	int i, j, k;
	double res;
	int half_pulse_num, dim2;
	double tau, two_tau;
	CLD res_amplitude;

	CLD ** exp_H1_tau = new CLD * [blk_num];
	CLD ** exp_H1_2tau = new CLD * [blk_num];
	CLD ** exp_H2_tau = new CLD * [blk_num];
	CLD ** exp_H2_2tau = new CLD * [blk_num];
	CLD ** mat_D1 = new CLD * [blk_num];
	CLD ** mat_D2 = new CLD * [blk_num];
	CLD ** mat_D2_N = new CLD * [blk_num];
	CLD ** res_mat = new CLD * [blk_num];
	CLD ** res_mat1 = new CLD * [blk_num];

	for(i = 0; i < blk_num; i++)
	{
		dim2 = blk_dim[i] * blk_dim[i];
		exp_H1_tau[i] = new CLD [dim2];
		exp_H1_2tau[i] = new CLD [dim2];
		exp_H2_tau[i] = new CLD [dim2];
		exp_H2_2tau[i] = new CLD [dim2];
		mat_D1[i] = new CLD[dim2];
		mat_D2[i] = new CLD[dim2];
		mat_D2_N[i] = new CLD[dim2];
		res_mat[i] = new CLD[dim2];
		res_mat1[i] = new CLD[dim2];
	}

	int highest_pow = PULSE_SYMMETRY_MODE[0];
	CLD *** mat_D2_n_list = new CLD ** [ highest_pow ];

	for(i = 0; i < highest_pow; i++)
	{
		mat_D2_n_list[i] = new CLD * [blk_num];
		for(j = 0; j < blk_num; j ++)
		{
			dim2 = blk_dim[j] * blk_dim[j];
			mat_D2_n_list[i][j] = new CLD [dim2];
		}
	}

	half_pulse_num = pulse_num / 2;	// = N for 2N pulses
	tau = timeseq[0];
	two_tau = timeseq[1];

	for(i = 0; i < blk_num; i++)
	{
		MKL_i_t_Hermite_Matrix_Exp(blk_dim[i], eigen_val1[i], eigen_vec1[i], -tau, exp_H1_tau[i]);
		MKL_i_t_Hermite_Matrix_Exp(blk_dim[i], eigen_val2[i], eigen_vec2[i], -tau, exp_H2_tau[i]);

		MKL_mat_mult(blk_dim[i], "NoTranspose", "NoTranspose", COMPLEX_ONE, exp_H1_tau[i], exp_H1_tau[i], COMPLEX_ZERO, exp_H1_2tau[i]);
		MKL_mat_mult(blk_dim[i], "NoTranspose", "NoTranspose", COMPLEX_ONE, exp_H2_tau[i], exp_H2_tau[i], COMPLEX_ZERO, exp_H2_2tau[i]);

		MKL_mat_mult(blk_dim[i], "NoTranspose", "NoTranspose", COMPLEX_ONE, exp_H1_tau[i], exp_H2_tau[i], COMPLEX_ZERO, mat_D1[i]);
		MKL_mat_mult(blk_dim[i], "NoTranspose", "NoTranspose", COMPLEX_ONE, exp_H1_2tau[i], exp_H2_2tau[i], COMPLEX_ZERO, mat_D2[i]);
	}

	for(i = 0; i < blk_num; i++)
		MKL_matrix_copy('N', blk_dim[i], COMPLEX_ONE, mat_D2[i], mat_D2_n_list[0][i]);

	for(i = 0; i < highest_pow - 1; i ++)
	{
		for(j = 0; j < blk_num; j++)
		{
			MKL_mat_mult(blk_dim[j], "NoTranspose", "NoTranspose", COMPLEX_ONE,
					mat_D2_n_list[i][j], mat_D2_n_list[i][j], COMPLEX_ZERO, mat_D2_n_list[i+1][j]);
		}
	}

	for(i = 0; i < blk_num; i++)
		MKL_matrix_copy('N', blk_dim[i], COMPLEX_ONE, mat_D2_n_list[ highest_pow - 1 ][i], mat_D2_N[i]);

	for(i = 0; i < highest_pow - 1; i ++)
	{
		if(PULSE_SYMMETRY_MODE[i+1] == 1)
		{
			for(j = 0; j < blk_num; j++)
			{
				MKL_mat_mult(blk_dim[j], "NoTranspose", "NoTranspose", COMPLEX_ONE,
						mat_D2_n_list[i][j], mat_D2_N[j], COMPLEX_ZERO, res_mat1[j]);
				MKL_matrix_copy('N', blk_dim[j], COMPLEX_ONE, res_mat1[j], mat_D2_N[j]);
			}
		}
	}

	for(i = 0; i < blk_num; i++)
	{
		MKL_mat_mult(blk_dim[i], "ConjugateTranspose", "NoTranspose", COMPLEX_ONE,
				mat_D1[i], mat_D2_N[i], COMPLEX_ZERO, res_mat[i]);

		MKL_mat_mult(blk_dim[i], "ConjugateTranspose", "NoTranspose", COMPLEX_ONE,
				mat_D2_N[i], res_mat[i], COMPLEX_ZERO, res_mat1[i]);

		MKL_mat_mult(blk_dim[i], "NoTranspose", "NoTranspose", COMPLEX_ONE,
				mat_D1[i], res_mat1[i], COMPLEX_ZERO, res_mat[i]);
	}

	for(i = 0; i < highest_pow; i ++)
	{
		for(j = 0; j < blk_num;j++)
			delete [] mat_D2_n_list[i][j];
		delete[] mat_D2_n_list[i];
	}
	delete[] mat_D2_n_list;


	res_amplitude = COMPLEX_ZERO;
	for(i = 0; i < blk_num; i++)
	{
		for(j = 0; j < blk_dim[i]; j ++)
			res_amplitude += res_mat[i][j * blk_dim[i] + j];

		delete [] exp_H1_tau[i];
		delete [] exp_H1_2tau[i];
		delete [] exp_H2_tau[i];
		delete [] exp_H2_2tau[i];
		delete [] mat_D1[i];
		delete [] mat_D2[i];
		delete [] mat_D2_N[i];
		delete [] res_mat[i];
		delete [] res_mat1[i];
	}
	delete [] exp_H1_tau;
	delete [] exp_H1_2tau;
	delete [] exp_H2_tau;
	delete [] exp_H2_2tau;
	delete [] mat_D1;
	delete [] mat_D2;
	delete [] mat_D2_N;
	delete [] res_mat;
	delete [] res_mat1;

//	res = abs(res_amplitude) / dim;
	res = res_amplitude.real() / dim;
	return res;
}

double cCLUSTER::PDD_algorithm_odd(double * timeseq)
{
	int i, j, k;
	double res;
	int half_pulse_num, dim2;
	double tau, two_tau;
	CLD res_amplitude;

	CLD ** exp_H1_tau = new CLD * [blk_num];
	CLD ** exp_H1_2tau = new CLD * [blk_num];
	CLD ** exp_H2_tau = new CLD * [blk_num];
	CLD ** exp_H2_2tau = new CLD * [blk_num];
	CLD ** mat_D1 = new CLD * [blk_num];
	CLD ** mat_D2 = new CLD * [blk_num];
	CLD ** mat_D2_N = new CLD * [blk_num];
	CLD ** mat_D_N = new CLD * [blk_num];
	CLD ** res_mat = new CLD * [blk_num];
	CLD ** res_mat1 = new CLD * [blk_num];

	for(i = 0; i < blk_num; i++)
	{
		dim2 = blk_dim[i] * blk_dim[i];
		exp_H1_tau[i] = new CLD [dim2];
		exp_H1_2tau[i] = new CLD [dim2];
		exp_H2_tau[i] = new CLD [dim2];
		exp_H2_2tau[i] = new CLD [dim2];
		mat_D1[i] = new CLD[dim2];
		mat_D2[i] = new CLD[dim2];
		mat_D2_N[i] = new CLD[dim2];
		mat_D_N[i] = new CLD[dim2];
		res_mat[i] = new CLD[dim2];
		res_mat1[i] = new CLD[dim2];
	}

	int highest_pow = PULSE_SYMMETRY_MODE[0];
	CLD *** mat_D2_n_list = new CLD ** [ highest_pow ];

	for(i = 0; i < highest_pow; i++)
	{
		mat_D2_n_list[i] = new CLD * [blk_num];
		for(j = 0; j < blk_num; j ++)
		{
			dim2 = blk_dim[j] * blk_dim[j];
			mat_D2_n_list[i][j] = new CLD [dim2];
		}
	}

	half_pulse_num = pulse_num / 2 + 1;	// = N for 2N-1 pulses
	tau = timeseq[0];
	two_tau = timeseq[1];

	for(i = 0; i < blk_num; i++)
	{
		MKL_i_t_Hermite_Matrix_Exp(blk_dim[i], eigen_val1[i], eigen_vec1[i], -tau, exp_H1_tau[i]);
		MKL_i_t_Hermite_Matrix_Exp(blk_dim[i], eigen_val2[i], eigen_vec2[i], -tau, exp_H2_tau[i]);

		MKL_mat_mult(blk_dim[i], "NoTranspose", "NoTranspose", COMPLEX_ONE, exp_H1_tau[i], exp_H1_tau[i], COMPLEX_ZERO, exp_H1_2tau[i]);
		MKL_mat_mult(blk_dim[i], "NoTranspose", "NoTranspose", COMPLEX_ONE, exp_H2_tau[i], exp_H2_tau[i], COMPLEX_ZERO, exp_H2_2tau[i]);

		MKL_mat_mult(blk_dim[i], "NoTranspose", "NoTranspose", COMPLEX_ONE, exp_H1_tau[i], exp_H2_tau[i], COMPLEX_ZERO, mat_D1[i]);
		MKL_mat_mult(blk_dim[i], "NoTranspose", "NoTranspose", COMPLEX_ONE, exp_H1_2tau[i], exp_H2_2tau[i], COMPLEX_ZERO, mat_D2[i]);
	}

	for(i = 0; i < blk_num; i++)
		MKL_matrix_copy('N', blk_dim[i], COMPLEX_ONE, mat_D2[i], mat_D2_n_list[0][i]);

	for(i = 0; i < highest_pow - 1; i ++)
	{
		for(j = 0; j < blk_num; j++)
		{
			MKL_mat_mult(blk_dim[j], "NoTranspose", "NoTranspose", COMPLEX_ONE,
					mat_D2_n_list[i][j], mat_D2_n_list[i][j], COMPLEX_ZERO, mat_D2_n_list[i+1][j]);
		}
	}

	for(i = 0; i < blk_num; i++)
		MKL_matrix_copy('N', blk_dim[i], COMPLEX_ONE, mat_D2_n_list[ highest_pow - 1 ][i], mat_D2_N[i]);

	for(i = 0; i < highest_pow - 1; i ++)
	{
		if(PULSE_SYMMETRY_MODE[i+1] == 1)
		{
			for(j = 0; j < blk_num; j++)
			{
				MKL_mat_mult(blk_dim[j], "NoTranspose", "NoTranspose", COMPLEX_ONE,
						mat_D2_n_list[i][j], mat_D2_N[j], COMPLEX_ZERO, res_mat1[j]);

				MKL_matrix_copy('N', blk_dim[j], COMPLEX_ONE, res_mat1[j], mat_D2_N[j]);
			}
		}
	}

	for(i = 0; i < blk_num; i++)
	{
		MKL_mat_mult(blk_dim[i], "NoTranspose", "NoTranspose", COMPLEX_ONE,
				mat_D1[i], mat_D2_N[i], COMPLEX_ZERO, res_mat1[i]);

		MKL_mat_mult(blk_dim[i], "ConjugateTranspose", "NoTranspose", COMPLEX_ONE,
				mat_D2[i], res_mat1[i], COMPLEX_ZERO, res_mat[i]);

		MKL_mat_mult(blk_dim[i], "ConjugateTranspose", "NoTranspose", COMPLEX_ONE,
				mat_D2_N[i], res_mat[i], COMPLEX_ZERO, res_mat1[i]);

		MKL_mat_mult(blk_dim[i], "NoTranspose", "NoTranspose", COMPLEX_ONE,
				mat_D1[i], res_mat1[i], COMPLEX_ZERO, res_mat[i]);
	}

	for(i = 0; i < highest_pow; i ++)
	{
		for(j = 0; j < blk_num;j++)
			delete [] mat_D2_n_list[i][j];
		delete[] mat_D2_n_list[i];
	}
	delete[] mat_D2_n_list;


	res_amplitude = COMPLEX_ZERO;
	for(i = 0; i < blk_num; i++)
	{
		for(j = 0; j < blk_dim[i]; j ++)
			res_amplitude += res_mat[i][j * blk_dim[i] + j];

		delete [] exp_H1_tau[i];
		delete [] exp_H1_2tau[i];
		delete [] exp_H2_tau[i];
		delete [] exp_H2_2tau[i];
		delete [] mat_D1[i];
		delete [] mat_D2[i];
		delete [] mat_D2_N[i];
		delete [] mat_D_N[i];
		delete [] res_mat[i];
		delete [] res_mat1[i];
	}
	delete [] exp_H1_tau;
	delete [] exp_H1_2tau;
	delete [] exp_H2_tau;
	delete [] exp_H2_2tau;
	delete [] mat_D1;
	delete [] mat_D2;
	delete [] mat_D2_N;
	delete [] mat_D_N;
	delete [] res_mat;
	delete [] res_mat1;

//	res = abs(res_amplitude) / dim;
	res = res_amplitude.real() / dim;
	return res;
}

double cCLUSTER::Generic_algorithm(double * timeseq)
{
	int i, j, k;
	int interval_num, count, count1;
	double sign, e_j, res = 0.0;
	CLD res_amplitude;


	//non-symmetric pulse
	interval_num = pulse_num + 1 ;
	count = 0;	count1 = 0;
	for (i = 0; i < interval_num; i ++)
	{
		count1 >= interval_num ? sign = 1.: sign = -1.;
		for(k = 0; k < blk_num; k ++)
		{
			for(j = 0; j < blk_dim[k]; j ++)
			{
				i_phase_fact[count] = CLD(0., sign * timeseq[2 * i] * eigen_val1[k][j]);
				count ++;
			}
		}
		count1 ++;

		count1 >= interval_num ? sign = 1.: sign = -1.;
		for(k = 0; k < blk_num; k ++)
		{
			for(j = 0; j < blk_dim[k]; j ++)
			{
				i_phase_fact[count] = CLD(0., sign * timeseq[2 * i + 1] * eigen_val2[k][j]);
				count ++;
			}
		}
		count1 ++;
	}
	vzExp(count,i_phase_fact,exp_eigen_val);

//		for(i = 0; i < count; i ++)
//		{
//			cout << i << ": " << i_phase_fact[i] << "; " << exp_eigen_val[i] << endl;
//		}
//		cout << endl;

	for(i = 0; i < blk_num; i ++)
	{
		MKL_matrix_copy('N', blk_dim[i], COMPLEX_ONE, right_mat0[i], right_mat[i]);
		MKL_matrix_copy('N', blk_dim[i], COMPLEX_ONE, left_mat0[i], left_mat[i]);
	}
	for(i = 0; i < pulse_num; i ++)
	{
		row_transform(exp_eigen_val + dim * 2 * i);
		W_Mat_transform(1);
		row_transform(exp_eigen_val + dim * (2 * i + 1));
		W_Mat_transform(2);

	}

	row_transform(exp_eigen_val + dim * 2 * pulse_num);
	W_Mat_transform(1);
	row_transform(exp_eigen_val + dim * (2 * pulse_num + 1));

	int dim2;
	CLD ** res_mat = new CLD *[blk_num];

	res_amplitude = COMPLEX_ZERO;
	for(i = 0; i < blk_num; i++)
	{
		dim2 = blk_dim[i] * blk_dim[i];
		res_mat[i] = new CLD [dim2];
		MKL_mat_mult(blk_dim[i], "NoTranspose", "NoTranspose", COMPLEX_ONE, left_mat0[i], right_mat[i], COMPLEX_ZERO, res_mat[i]);
		for(j = 0; j < blk_dim[i]; j ++)
			res_amplitude += res_mat[i][j * blk_dim[i] + j];
	}
	for(i = 0; i < blk_num; i++)
		delete [] res_mat[i];
	delete [] res_mat;

//	res = abs(res_amplitude) / dim;
	res = res_amplitude.real() / dim;
	return res;
}

void cCLUSTER::row_transform(CLD * exp_vals)
{

	int i;
	for(i = 0; i < blk_num; i ++)
	{
		MKL_matrix_complex_row_transform(blk_dim[i], exp_vals, right_mat[i]);
		exp_vals += blk_dim[i];
	}

}

void cCLUSTER::W_Mat_transform(int isTrans)
{
	int i, dim2;
	CLD ** res = new CLD * [blk_num];

    char * transa, * transb;
	isTrans == 1 ? transa = "NoTranspose" : transa = "ConjugateTranspose";
    transb = "NoTranspose";

	for(i = 0; i < blk_num; i++)
	{
		dim2 = blk_dim[i] * blk_dim[i];
		res[i] = new CLD [dim2];
		MKL_mat_mult(blk_dim[i], transa, transb, COMPLEX_ONE, W_matrix[i], right_mat[i], COMPLEX_ZERO, res[i]);

		MKL_matrix_copy('N', blk_dim[i], COMPLEX_ONE, res[i], right_mat[i]);
	}

//	for(i = 0; i < blk_num; i++)
//		cout_C_Hami(right_mat, blk_num, blk_dim, 1, 1);

	for(i = 0; i < blk_num; i++)
		delete [] res[i];
	delete [] res;

}

void cCLUSTER::set_localfield()
{
	int i, idx_i;
	double z;

	for (i = 0; i < nspin; i++)
	{
		idx_i = spin_index[i];

		z	= spinlist[ idx_i ].Zeeman();

		if(IsDriving)
		{
			locX  [i] = DrvRabi;
			locY  [i] = 0.;
			locZ  [i] = DrvDetuning;
		}
		else
		{
			locX  [i] = 0.;
			locY  [i] = 0.;
			locZ  [i] = z;
		}

	}
}
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//Correlation function

void cCLUSTER::initialize_Correlation(cSPIN * spins, int num_spin, int * spinIdx, int npulse, int total_idx, clusterHamiltonian hamiform)
{
	int i, j, k, idx_i, idx_j, idx_ij;
	assert(num_spin <= MAX_SPIN_IN_CLST);

	spinlist		= spins;
	nspin			= num_spin;
	spin_index		= spinIdx;
	pulse_num		= npulse;
	clst_idx		= total_idx;

	sINT_TENSOR1 dip;

	//	dim = 1;
	dim = hamiform.FullDim();
	for (i = 0; i < nspin; i++)
	{
//		dim *= 2;
		idx_i = spin_index[i];

		hf1x[i] = spinlist[idx_i].hf(1).x;
		hf1y[i] = spinlist[idx_i].hf(1).y;
		hf1z[i] = spinlist[idx_i].hf(1).z;

		hf2x[i] = spinlist[idx_i].hf(2).x;
		hf2y[i] = spinlist[idx_i].hf(2).y;
		hf2z[i] = spinlist[idx_i].hf(2).z;

		for (j = 0; j < nspin; j++)
		{
			idx_j = spin_index[j];
			idx_ij = idx_i * NBATH + idx_j;

			dip = calc_dipole_int(idx_i, idx_j);

			if(IsMediated)
			{
			xx1 [i * nspin + j] = dip.xx + MEDIATED_INT1[idx_ij + 0];
			yy1 [i * nspin + j] = dip.yy + MEDIATED_INT1[idx_ij + 1];
			zz1 [i * nspin + j] = dip.zz + MEDIATED_INT1[idx_ij + 2];

			xy1 [i * nspin + j] = dip.xy + MEDIATED_INT1[idx_ij + 3];
			yz1 [i * nspin + j] = dip.yz + MEDIATED_INT1[idx_ij + 4];
			zx1 [i * nspin + j] = dip.zx + MEDIATED_INT1[idx_ij + 5];

			yx1 [i * nspin + j] = dip.yx + MEDIATED_INT1[idx_ij + 6];
			zy1 [i * nspin + j] = dip.zy + MEDIATED_INT1[idx_ij + 7];
			xz1 [i * nspin + j] = dip.xz + MEDIATED_INT1[idx_ij + 8];

			// xx1==xx2...
			xx2 [i * nspin + j] = dip.xx + MEDIATED_INT2[idx_ij + 0];
			yy2 [i * nspin + j] = dip.yy + MEDIATED_INT2[idx_ij + 1];
			zz2 [i * nspin + j] = dip.zz + MEDIATED_INT2[idx_ij + 2];

			xy2 [i * nspin + j] = dip.xy + MEDIATED_INT2[idx_ij + 3];
			yz2 [i * nspin + j] = dip.yz + MEDIATED_INT2[idx_ij + 4];
			zx2 [i * nspin + j] = dip.zx + MEDIATED_INT2[idx_ij + 5];

			yx2 [i * nspin + j] = dip.yx + MEDIATED_INT2[idx_ij + 6];
			zy2 [i * nspin + j] = dip.zy + MEDIATED_INT2[idx_ij + 7];
			xz2 [i * nspin + j] = dip.xz + MEDIATED_INT2[idx_ij + 8];
			}
			else
			{
				xx1 [i * nspin + j] = dip.xx;
				yy1 [i * nspin + j] = dip.yy;
				zz1 [i * nspin + j] = dip.zz;

				xy1 [i * nspin + j] = dip.xy;
				yz1 [i * nspin + j] = dip.yz;
				zx1 [i * nspin + j] = dip.zx;

				yx1 [i * nspin + j] = dip.yx;
				zy1 [i * nspin + j] = dip.zy;
				xz1 [i * nspin + j] = dip.xz;

				// xx1==xx2...
				xx2 [i * nspin + j] = dip.xx;
				yy2 [i * nspin + j] = dip.yy;
				zz2 [i * nspin + j] = dip.zz;

				xy2 [i * nspin + j] = dip.xy;
				yz2 [i * nspin + j] = dip.yz;
				zx2 [i * nspin + j] = dip.zx;

				yx2 [i * nspin + j] = dip.yx;
				zy2 [i * nspin + j] = dip.zy;
				xz2 [i * nspin + j] = dip.xz;
			}
		}
	}
	ZEROS[0] = 0.0;

	set_localfield();

	hm          = hamiform;
	blk_num     = hm.BlockNumer();
	blk_dim     = hm.HamiDim();

	C_HAMI_0 = new CLD * [blk_num];     eigen_val0 = new double * [blk_num];    eigen_vec0 = new CLD * [blk_num];
	C_HAMI_1 = new CLD * [blk_num];     eigen_val1 = new double * [blk_num];    eigen_vec1 = new CLD * [blk_num];
	C_HAMI_2 = new CLD * [blk_num];     eigen_val2 = new double * [blk_num];    eigen_vec2 = new CLD * [blk_num];
	H_B_Mat  = new CLD * [blk_num];		eigen_val_H_B = new double * [blk_num];    eigen_vec_H_B = new CLD * [blk_num];
	b_Mat	 = new CLD * [blk_num];		eigen_val_b = new double * [blk_num];    eigen_vec_b = new CLD * [blk_num];
	W_matrix = new CLD * [blk_num];

	left_mat0  = new CLD * [blk_num];
	left_mat   = new CLD * [blk_num];
	right_mat0 = new CLD * [blk_num];
	right_mat  = new CLD * [blk_num];

	int n_dim1, n_dim2, n_dim = 0;
	for(i = 0; i < blk_num; i++)
	{
		n_dim += blk_dim[i];
		n_dim1 = blk_dim[i];
		n_dim2 = n_dim1 * n_dim1;
		C_HAMI_0[i] = new CLD [ n_dim2 ];   eigen_val0[i] = new double [n_dim1];    eigen_vec0[i] = new CLD [n_dim2];
		C_HAMI_1[i] = new CLD [ n_dim2 ];   eigen_val1[i] = new double [n_dim1];    eigen_vec1[i] = new CLD [n_dim2];
		C_HAMI_2[i] = new CLD [ n_dim2 ];   eigen_val2[i] = new double [n_dim1];    eigen_vec2[i] = new CLD [n_dim2];
		H_B_Mat[i]  = new CLD [ n_dim2 ];	eigen_val_H_B[i] =  new double [n_dim1];eigen_vec_H_B[i] = new CLD [n_dim2];
		b_Mat[i]  = new CLD [ n_dim2 ];		eigen_val_b[i] =  new double [n_dim1];  eigen_vec_b[i] = new CLD [n_dim2];
		W_matrix[i] = new CLD [ n_dim2 ];

		left_mat0[i]  = new CLD [n_dim2];
		right_mat0[i] = new CLD [n_dim2];
		left_mat[i]   = new CLD [n_dim2];
		right_mat[i]  = new CLD [n_dim2];
	}
	assert(n_dim == dim);

	i_phase_fact = new CLD [ dim * 2 * (pulse_num + 1)];
	exp_eigen_val = new CLD [dim * 2 * (pulse_num + 1)];
}

void cCLUSTER::calc_Hamiltonian_Correlation(int sign)
{
	int i, j;

        CLD     ** C_HAMI;
        double  ** eigenval;
	CLD 	** eigenvec;

	double * xx, * yy, * zz, * xy, * yz, * zx, * yx, * zy, * xz;

	if(sign == 1)
	{
		for (i = 0; i < nspin; i++)
		{
			Bx[i] = locX[i] + hf1x[i];
			By[i] = locY[i] + hf1y[i];
			Bz[i] = locZ[i] + hf1z[i];

		}
		xx = xx1; yy = yy1; zz = zz1;
		xy = xy1; yz = yz1; zx = zx1;
		yx = yx1; zy = zy1; xz = xz1;

		C_HAMI = C_HAMI_1;
		eigenval = eigen_val1;
		eigenvec = eigen_vec1;
	}
	else if(sign == 2)
	{
		for (i = 0; i < nspin; i++)
		{
			Bx[i] = locX[i] + hf2x[i];
			By[i] = locY[i] + hf2y[i];
			Bz[i] = locZ[i] + hf2z[i];
		}
		xx = xx2; yy = yy2; zz = zz2;
		xy = xy2; yz = yz2; zx = zx2;
		yx = yx2; zy = zy2; xz = xz2;

		C_HAMI = C_HAMI_2;
		eigenval = eigen_val2;
		eigenvec = eigen_vec2;

	}
	else // sign==3
	{
		for (i = 0; i < nspin; i++)
		{
			Bx[i] = locX[i] + 0.5 * (hf1x[i] + hf2x[i]);
			By[i] = locY[i] + 0.5 * (hf1y[i] + hf2y[i]);
			Bz[i] = locZ[i] + 0.5 * (hf1z[i] + hf2z[i]);
		}
		xx = xx2; yy = yy2; zz = zz2;
		xy = xy2; yz = yz2; zx = zx2;
		yx = yx2; zy = zy2; xz = xz2;

		C_HAMI = C_HAMI_0;
		eigenval = eigen_val0;
		eigenvec = eigen_vec0;
	}

	INTERACTIONS[0] = ZEROS;
	INTERACTIONS[1] = Bx; INTERACTIONS[2] = By; INTERACTIONS[3] = Bz;
	INTERACTIONS[4] = xx; INTERACTIONS[5] = yy; INTERACTIONS[6] = zz;
	INTERACTIONS[7] = xy; INTERACTIONS[8] = yz; INTERACTIONS[9] = zx;
	INTERACTIONS[10] = yx; INTERACTIONS[11] = zy; INTERACTIONS[12] = xz;

	double re, im;
	for(i = 0; i < blk_num; i++)
	{
		for(j = 0; j < blk_dim[i]*blk_dim[i]; j++)
		{
			re = hm.ElementVal_R(INTERACTIONS, i, j);
			im = hm.ElementVal_I(INTERACTIONS, i, j);
			C_HAMI[i][j] = CLD(re, im);
		}
	}

//	if(sign == 1)
//		diagnosis_CLD("tst1.dat", C_HAMI[0],  blk_dim[0],  blk_dim[0]);
//	else
//		diagnosis_CLD("tst2.dat", C_HAMI[0],  blk_dim[0],  blk_dim[0]);

}

void cCLUSTER::pre_Signal_Correlation()
{
    int i, j, k, dim_i;
    for(i = 0; i < blk_num; i++)
    {
    	dim_i = blk_dim[i];
        for(j = 0; j < dim_i; j++)
        {
        	for(k = 0; k < dim_i; k++)
        	{
        		H_B_Mat[i][ j * dim_i + k] = 0.5 * (C_HAMI_1[i][ j * dim_i + k] + C_HAMI_2[i][ j * dim_i + k]);
        		b_Mat [i][ j * dim_i + k]  =
        				j <= k ?
        				(C_HAMI_1[i][ j * dim_i + k] - C_HAMI_2[i][ j * dim_i + k]):
        				conj(C_HAMI_1[i][ k * dim_i + j] - C_HAMI_2[i][ k * dim_i + j]);
        	}
        }
    }

    CLD     ** C_HAMI, ** eigenvec;
    double  ** eigenval;

//    diagnosis_CLD("HbMat.dat", H_B_Mat[0],  blk_dim[0],  blk_dim[0]);
//    diagnosis_CLD("bMat.dat", b_Mat[0],  blk_dim[0],  blk_dim[0]);

    C_HAMI = H_B_Mat;    eigenval = eigen_val_H_B;    eigenvec = eigen_vec_H_B;
    for(i = 0; i < blk_num; i++)
        MKL_hermite_mat_diag(blk_dim[i], C_HAMI[i], eigenval[i],eigenvec[i]);

//    diagnosis_CLD("EigenVec.dat", eigenvec[0],  blk_dim[0],  blk_dim[0]);

    CLD		** W_matrix_temp;
    W_matrix_temp =  new CLD * [blk_num];
    for(i = 0; i < blk_num; i++)
    	W_matrix_temp[i] = new CLD[ blk_dim[i] * blk_dim[i] ];

    char * transa, * transb;
    transa = "ConjugateTranspose";
    transb = "NoTranspose";

    for(i = 0; i < blk_num; i++)
    {
        MKL_mat_mult(blk_dim[i], transa, transb, COMPLEX_ONE, eigenvec[i], b_Mat[i], COMPLEX_ZERO, W_matrix_temp[i]);
        MKL_mat_mult(blk_dim[i], transb, transb, COMPLEX_ONE, W_matrix_temp[i], eigenvec[i], COMPLEX_ZERO, W_matrix[i]);
    }

    for(i = 0; i < blk_num; i++)
    	delete [] W_matrix_temp[i];
    delete [] W_matrix_temp; W_matrix_temp =NULL;

//    diagnosis_CLD("Wmat.dat", W_matrix[0],  blk_dim[0],  blk_dim[0]);
}

double cCLUSTER::CorrelationFunc(double tval)
{
	int i, j, k, dim_i;
	double res, weight, w_j, w_k;

	res = 0.0;
	for(i = 0; i < blk_num; i++)
	{
		dim_i = blk_dim[i];
		for(j = 0; j < dim_i; j++)
		{
			res += (W_matrix[i][ j*dim_i + j ] * W_matrix[i][ j*dim_i + j ]).real();
			for(k = 0; k < j; k++)
			{
				weight = (W_matrix[i][ j*dim_i + k ] * W_matrix[i][ k*dim_i + j ]).real();
				w_j = eigen_val_H_B[i][j];	w_k = eigen_val_H_B[i][k];
				res += 2.0 * weight * cos( tval* (w_j - w_k) );
			}
		}
	}
	res /= dim;

	return res;
}

double cCLUSTER::NoiseSpectrum(double omega)
{
	int i, j, k, dim_i;
	double res, weight0, weight, w_j, w_k;

	res = 0.0;
	for(i = 0; i < blk_num; i++)
	{
		dim_i = blk_dim[i];
		for(j = 0; j < dim_i; j++)
		{
//			weight0 = (W_matrix[i][ j*dim_i + j ] * W_matrix[i][ j*dim_i + j ]).real();
//			res += weight0 * Lorentzian(omega, 0.0, EXP_DECAY_RATE);
			for(k = 0; k < j; k++)
			{
				weight = (W_matrix[i][ j*dim_i + k ] * W_matrix[i][ k*dim_i + j ]).real();
				w_j = eigen_val_H_B[i][j];	w_k = eigen_val_H_B[i][k];
//				if(fabs(w_j - w_k) > 0.5 * EXP_DECAY_RATE)
					res += 2.0 * weight *  Lorentzian(omega, w_j - w_k, EXP_DECAY_RATE);
			}
		}
	}
	res /= dim;

	return res;
}

double cCLUSTER::GaussianEcho(double tval)
{
	int i, j, k, dim_i;
	double res, int_res, weight, w_j, w_k, dw;

	int_res = 0.0;
	for(i = 0; i < blk_num; i++)
	{
		dim_i = blk_dim[i];
		for(j = 0; j < dim_i; j++)
		{
//			int_res += (W_matrix[i][ j*dim_i + j ] * W_matrix[i][ j*dim_i + j ]).real();
			for(k = 0; k < j; k++)
			{
				weight = (W_matrix[i][ j*dim_i + k ] * W_matrix[i][ k*dim_i + j ]).real();

				w_j = eigen_val_H_B[i][j];	w_k = eigen_val_H_B[i][k];
				dw = w_j - w_k;

				int_res += 2.0 * weight * int_cos_function(dw, tval);
			}
		}
	}
	int_res /= dim;
	res = exp(-0.5 * int_res);

	return res;
//
//	int i, j;
//	double res, weight, w_i, w_j;
//
//	res = 0;
//	for (i = 0; i < dim ; i++)
//	{
//		w_i = gsl_vector_get(eigen_val0, i);
//		weight = gsl_complex_abs2( gsl_matrix_complex_get(WEIGHT_MAT, i, i) );
//		res += weight * int_cos_func(0.0, tval, pulse_num);
//		for ( j = 0; j < i; j++)
//		{
//			w_j = gsl_vector_get(eigen_val0, j);
//
//			weight = gsl_complex_abs2( gsl_matrix_complex_get(WEIGHT_MAT, i, j) );
//			res += 2.0 * weight * int_cos_func(w_i - w_j, tval, pulse_num);
//		}
//	}
//
//	res /= dim;
//	return res;
}

double cCLUSTER::int_cos_func(double delta_w, double tval, int npulse)
{
    double fact, delta_t, detla_w2, res;

    int i, j;
    double sign_i, sign_j;

    res = 0.0;
    detla_w2 = delta_w*delta_w;

//    if(fabs(delta_w) < 1.0E-6 && npulse != 0) return res;
//    if(fabs(delta_w) < 1.0E-6 && npulse == 0) return tval*tval;

    sign_i = 1.0;
    for(i = 0; i < npulse + 2; i++)
    {
        sign_j = 1.0;
        for(j = 0; j < i; j++)
        {
            if( (j == 0 && i != npulse + 1) || ( j != 0 && i == npulse + 1)  )
            {
                fact = 4.0;
            }
            else if( j == 0 && i == npulse + 1)
            {
                fact = 2.0;
            }
            else
            {
                fact = 8.0;
            }

            delta_t = tval* ( PULSE_TIMING[i] - PULSE_TIMING[j] );

            if( fabs (tval*delta_w) > 1.0E-5)
                res += sign_i * sign_j * fact * cos(delta_t * delta_w)/detla_w2;
            else
                res += sign_i * sign_j * fact * (-0.5*delta_t*delta_t);

            sign_j *= -1.0;
        }

        if(fabs(tval*delta_w) > 1.0E-5)
        {
            if(i == 0 || i == npulse + 1)
                res += 1.0/detla_w2;
            else
                res += 4.0/detla_w2;
        }

        sign_i *= -1.0;
    }


    return res;
}

double cCLUSTER::int_cos_function(double delta_w, double tval)
{
    double fact, delta_t, delta_w2, res;

    int i, j;
    double sign_i, sign_j;

    res = 0.0;
    delta_w2 = delta_w*delta_w;

    sign_i = 1.0;
    for(i = 0; i < PULSE_NUM + 2; i++)
    {
        sign_j = 1.0;
        for(j = 0; j < PULSE_NUM + 2; j++)
        {
        	if( 	(i == 0 && j == 0)
        		 || (i == 0 && j == PULSE_NUM + 1)
        		 || (i == PULSE_NUM + 1 && j == PULSE_NUM + 1)
        		 || (i == PULSE_NUM + 1 && j == 0) )// four corners
        	{
        		fact = 1.0;
        	}
        	else if( 	(j == 0 && i != PULSE_NUM + 1 && i != 0)
        			 || (j == PULSE_NUM +1 && i != PULSE_NUM + 1 && i != 0)
        			 || (i == 0 && j != PULSE_NUM + 1 && j != 0)
        			 || (i == PULSE_NUM +1 && j != PULSE_NUM + 1 && j != 0) )
            {
                fact = 2.0;
            }
            else
            {
                fact = 4.0;
            }

            delta_t = tval* ( PULSE_TIMING[i] - PULSE_TIMING[j] );

            res += sign_i * sign_j * fact * cos(delta_t * delta_w);

            sign_j *= -1.0;
        }
        sign_i *= -1.0;
    }

    res /= delta_w2;
    return res;
}

void cCLUSTER::clear_Correlation()
{
    int i;
    for(i = 0; i < blk_num; i++)
    {
        delete [] C_HAMI_0[i];  delete [] eigen_val0[i];    delete [] eigen_vec0[i];
        delete [] C_HAMI_1[i];  delete [] eigen_val1[i];    delete [] eigen_vec1[i];
        delete [] C_HAMI_2[i];  delete [] eigen_val2[i];    delete [] eigen_vec2[i];
        delete [] H_B_Mat[i]; 	delete [] eigen_val_H_B[i];	delete [] eigen_vec_H_B[i];
        delete [] b_Mat[i];		delete [] eigen_val_b[i];	delete [] eigen_vec_b[i];
        delete [] W_matrix[i];
        delete [] left_mat0[i];		delete [] left_mat[i];
        delete [] right_mat0[i];	delete [] right_mat[i];
    }
    delete [] C_HAMI_0; C_HAMI_0 = NULL;
    delete [] C_HAMI_1; C_HAMI_1 = NULL;
    delete [] C_HAMI_2; C_HAMI_2 = NULL;
    delete [] H_B_Mat; 	H_B_Mat = NULL;
    delete [] b_Mat; 	b_Mat = NULL;

    delete [] eigen_val0; eigen_val0 = NULL;
    delete [] eigen_val1; eigen_val1 = NULL;
    delete [] eigen_val2; eigen_val2 = NULL;
    delete [] eigen_val_H_B; eigen_val_H_B = NULL;
    delete [] eigen_val_b; eigen_val_b = NULL;

    delete [] eigen_vec0; eigen_vec0 = NULL;
    delete [] eigen_vec1; eigen_vec1 = NULL;
    delete [] eigen_vec2; eigen_vec2 = NULL;
    delete [] eigen_vec_H_B; eigen_vec_H_B = NULL;
    delete [] eigen_vec_b; eigen_vec_b = NULL;

    delete [] W_matrix; W_matrix =NULL;
    delete [] left_mat0; left_mat0 = NULL;
    delete [] right_mat0; right_mat0 = NULL;
    delete [] left_mat; left_mat = NULL;
    delete [] right_mat; right_mat = NULL;

    delete [] i_phase_fact; i_phase_fact = NULL;
    delete [] exp_eigen_val; exp_eigen_val = NULL;
}

////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//Single Sample

void cCLUSTER::initialize_SGL(cSPIN * spins, int num_spin, int * spinIdx, int npulse, int total_idx, clusterHamiltonian hamiform,  int clst_state_j)
{
	int i, j, k, idx_i, idx_j, idx_ij;
	assert(num_spin <= MAX_SPIN_IN_CLST);

	spinlist		= spins;
	nspin			= num_spin;
	spin_index		= spinIdx;
	pulse_num		= npulse;
	clst_idx		= total_idx;

	clster_state_J  = clst_state_j;

	sINT_TENSOR1 dip;

	double s;
	s = (IsMediated == 1 ? 1.0 : 0.0);

//	dim = 1;
	dim = hamiform.FullDim();
	for (i = 0; i < nspin; i++)
	{
//		dim *= 2;
		idx_i = spin_index[i];

		hf1x[i] = spinlist[idx_i].hf(1).x;
		hf1y[i] = spinlist[idx_i].hf(1).y;
		hf1z[i] = spinlist[idx_i].hf(1).z;

		hf2x[i] = spinlist[idx_i].hf(2).x;
		hf2y[i] = spinlist[idx_i].hf(2).y;
		hf2z[i] = spinlist[idx_i].hf(2).z;

		for (j = 0; j < nspin; j++)
		{
			idx_j = spin_index[j];
			idx_ij = idx_i * NBATH + idx_j;

			dip = calc_dipole_int(idx_i, idx_j);

			if(IsMediated)
			{
				xx1 [i * nspin + j] = dip.xx + MEDIATED_INT1[idx_ij + 0];
				yy1 [i * nspin + j] = dip.yy + MEDIATED_INT1[idx_ij + 1];
				zz1 [i * nspin + j] = dip.zz + MEDIATED_INT1[idx_ij + 2];

				xy1 [i * nspin + j] = dip.xy + MEDIATED_INT1[idx_ij + 3];
				yz1 [i * nspin + j] = dip.yz + MEDIATED_INT1[idx_ij + 4];
				zx1 [i * nspin + j] = dip.zx + MEDIATED_INT1[idx_ij + 5];

				yx1 [i * nspin + j] = dip.yx + MEDIATED_INT1[idx_ij + 6];
				zy1 [i * nspin + j] = dip.zy + MEDIATED_INT1[idx_ij + 7];
				xz1 [i * nspin + j] = dip.xz + MEDIATED_INT1[idx_ij + 8];

				// xx1==xx2...
				xx2 [i * nspin + j] = dip.xx + MEDIATED_INT2[idx_ij + 0];
				yy2 [i * nspin + j] = dip.yy + MEDIATED_INT2[idx_ij + 1];
				zz2 [i * nspin + j] = dip.zz + MEDIATED_INT2[idx_ij + 2];

				xy2 [i * nspin + j] = dip.xy + MEDIATED_INT2[idx_ij + 3];
				yz2 [i * nspin + j] = dip.yz + MEDIATED_INT2[idx_ij + 4];
				zx2 [i * nspin + j] = dip.zx + MEDIATED_INT2[idx_ij + 5];

				yx2 [i * nspin + j] = dip.yx + MEDIATED_INT2[idx_ij + 6];
				zy2 [i * nspin + j] = dip.zy + MEDIATED_INT2[idx_ij + 7];
				xz2 [i * nspin + j] = dip.xz + MEDIATED_INT2[idx_ij + 8];
			}
			else
			{
				xx1 [i * nspin + j] = dip.xx;
				yy1 [i * nspin + j] = dip.yy;
				zz1 [i * nspin + j] = dip.zz;

				xy1 [i * nspin + j] = dip.xy;
				yz1 [i * nspin + j] = dip.yz;
				zx1 [i * nspin + j] = dip.zx;

				yx1 [i * nspin + j] = dip.yx;
				zy1 [i * nspin + j] = dip.zy;
				xz1 [i * nspin + j] = dip.xz;

				// xx1==xx2...
				xx2 [i * nspin + j] = dip.xx;
				yy2 [i * nspin + j] = dip.yy;
				zz2 [i * nspin + j] = dip.zz;

				xy2 [i * nspin + j] = dip.xy;
				yz2 [i * nspin + j] = dip.yz;
				zx2 [i * nspin + j] = dip.zx;

				yx2 [i * nspin + j] = dip.yx;
				zy2 [i * nspin + j] = dip.zy;
				xz2 [i * nspin + j] = dip.xz;
			}
		}
	}
	ZEROS[0] = 0.0;

	set_localfield_SGL();

	hm          = hamiform;
	blk_num     = hm.BlockNumer();
	blk_dim     = hm.HamiDim();

	need_calc	= new int [blk_num];

	C_HAMI_0 = new CLD * [blk_num];     eigen_val0 = new double * [blk_num];    eigen_vec0 = new CLD * [blk_num];
	C_HAMI_1 = new CLD * [blk_num];     eigen_val1 = new double * [blk_num];    eigen_vec1 = new CLD * [blk_num];
	C_HAMI_2 = new CLD * [blk_num];     eigen_val2 = new double * [blk_num];    eigen_vec2 = new CLD * [blk_num];
	W_matrix = new CLD * [blk_num];

	left_vec0  = new CLD * [blk_num];
	left_vec   = new CLD * [blk_num];
	right_vec0 = new CLD * [blk_num];
	right_vec  = new CLD * [blk_num];

	int n_dim1, n_dim2, n_dim = 0;
	for(i = 0; i < blk_num; i++)
	{
		n_dim += blk_dim[i];
		n_dim1 = blk_dim[i];

		if( (clster_state_J >= n_dim - n_dim1) && (clster_state_J < n_dim) )
			need_calc[i] = clster_state_J - (n_dim - n_dim1);
		else
			need_calc[i] = -1;

		n_dim2 = n_dim1 * n_dim1;
		C_HAMI_0[i] = new CLD [ n_dim2 ];   eigen_val0[i] = new double [n_dim1];    eigen_vec0[i] = new CLD [n_dim2];
		C_HAMI_1[i] = new CLD [ n_dim2 ];   eigen_val1[i] = new double [n_dim1];    eigen_vec1[i] = new CLD [n_dim2];
		C_HAMI_2[i] = new CLD [ n_dim2 ];   eigen_val2[i] = new double [n_dim1];    eigen_vec2[i] = new CLD [n_dim2];
		W_matrix[i] = new CLD [ n_dim2 ];

		left_vec0[i]  = new CLD [n_dim1];
		right_vec0[i] = new CLD [n_dim1];
		left_vec[i]   = new CLD [n_dim1];
		right_vec[i]  = new CLD [n_dim1];
	}
	assert(n_dim == dim);

	i_phase_fact = new CLD [ dim * 2 * (pulse_num + 1)];
	exp_eigen_val = new CLD [dim * 2 * (pulse_num + 1)];

}

void cCLUSTER::clear_SGL()
{
    int i;
    for(i = 0; i < blk_num; i++)
    {
        delete [] C_HAMI_0[i];  delete [] eigen_val0[i];    delete [] eigen_vec0[i];
        delete [] C_HAMI_1[i];  delete [] eigen_val1[i];    delete [] eigen_vec1[i];
        delete [] C_HAMI_2[i];  delete [] eigen_val2[i];    delete [] eigen_vec2[i];
        delete [] W_matrix[i];
        delete [] left_vec0[i];		delete [] left_vec[i];
        delete [] right_vec0[i];	delete [] right_vec[i];
    }
    delete [] C_HAMI_0; C_HAMI_0 = NULL;
    delete [] C_HAMI_1; C_HAMI_1 = NULL;
    delete [] C_HAMI_2; C_HAMI_2 = NULL;

    delete [] eigen_val0; eigen_val0 = NULL;
    delete [] eigen_val1; eigen_val1 = NULL;
    delete [] eigen_val2; eigen_val2 = NULL;

    delete [] eigen_vec0; eigen_vec0 = NULL;
    delete [] eigen_vec1; eigen_vec1 = NULL;
    delete [] eigen_vec2; eigen_vec2 = NULL;

    delete [] W_matrix; W_matrix =NULL;
    delete [] left_vec0; left_vec0 = NULL;
    delete [] right_vec0; right_vec0 = NULL;
    delete [] left_vec; left_vec = NULL;
    delete [] right_vec; right_vec = NULL;

    delete [] i_phase_fact; i_phase_fact = NULL;
    delete [] exp_eigen_val; exp_eigen_val = NULL;
}

void cCLUSTER::calc_Hamiltonian_SGL(int sign)
{
	int i, j;

	CLD     ** C_HAMI;
	double  ** eigenval;
	CLD 	** eigenvec;

	double * xx, * yy, * zz, * xy, * yz, * zx, * yx, * zy, * xz;

	if(sign == 1)
	{
		for (i = 0; i < nspin; i++)
		{
			Bx[i] = locX[i] + hf1x[i];
			By[i] = locY[i] + hf1y[i];
			Bz[i] = locZ[i] + hf1z[i];

		}
		xx = xx1; yy = yy1; zz = zz1;
		xy = xy1; yz = yz1; zx = zx1;
		yx = yx1; zy = zy1; xz = xz1;

		C_HAMI = C_HAMI_1;
		eigenval = eigen_val1;
		eigenvec = eigen_vec1;
	}
	else if(sign == 2)
	{
		for (i = 0; i < nspin; i++)
		{
			Bx[i] = locX[i] + hf2x[i];
			By[i] = locY[i] + hf2y[i];
			Bz[i] = locZ[i] + hf2z[i];
		}
		xx = xx2; yy = yy2; zz = zz2;
		xy = xy2; yz = yz2; zx = zx2;
		yx = yx2; zy = zy2; xz = xz2;

		C_HAMI = C_HAMI_2;
		eigenval = eigen_val2;
		eigenvec = eigen_vec2;

	}
	else // sign==3
	{
		for (i = 0; i < nspin; i++)
		{
			Bx[i] = locX[i] + 0.5 * (hf1x[i] + hf2x[i]);
			By[i] = locY[i] + 0.5 * (hf1y[i] + hf2y[i]);
			Bz[i] = locZ[i] + 0.5 * (hf1z[i] + hf2z[i]);
		}
		xx = xx2; yy = yy2; zz = zz2;
		xy = xy2; yz = yz2; zx = zx2;
		yx = yx2; zy = zy2; xz = xz2;

		C_HAMI = C_HAMI_0;
		eigenval = eigen_val0;
		eigenvec = eigen_vec0;
	}

	INTERACTIONS[0] = ZEROS;
	INTERACTIONS[1] = Bx; INTERACTIONS[2] = By; INTERACTIONS[3] = Bz;
	INTERACTIONS[4] = xx; INTERACTIONS[5] = yy; INTERACTIONS[6] = zz;
	INTERACTIONS[7] = xy; INTERACTIONS[8] = yz; INTERACTIONS[9] = zx;
	INTERACTIONS[10] = yx; INTERACTIONS[11] = zy; INTERACTIONS[12] = xz;

	double re, im;
	for(i = 0; i < blk_num; i++)
	{
		if(need_calc[i] >= 0)
		{
			for(j = 0; j < blk_dim[i]*blk_dim[i]; j++)
			{
				re = hm.ElementVal_R(INTERACTIONS, i, j);
				im = hm.ElementVal_I(INTERACTIONS, i, j);
				C_HAMI[i][j] = CLD(re, im);
			}
		}
	}

//	cout_C_Hami(C_HAMI, blk_num, blk_dim, clst_idx, sign);

	for(i = 0; i < blk_num; i++)
	{
		if(need_calc[i] >= 0)
			MKL_hermite_mat_diag(blk_dim[i], C_HAMI[i], eigenval[i],eigenvec[i]);
	}
}

void cCLUSTER::pre_Signal_SGL()
{
    int i, j, dim_i;
    char * transa, * transb;

    transa = "ConjugateTranspose";
    transb = "NoTranspose";

    int state_id;
    for(i = 0; i < blk_num; i++)
    {
    	if(need_calc[i] >= 0)
    	{
    		state_id = need_calc[i];
    		dim_i = blk_dim[i];
    		MKL_mat_mult(dim_i, transa, transb, COMPLEX_ONE, eigen_vec2[i], eigen_vec1[i], COMPLEX_ZERO, W_matrix[i]);

    		for(j = 0; j < dim_i; j++)
    		{
    			left_vec[i][j]  = 		eigen_vec2[i][ dim_i * j + state_id];		// the state_id'th row of eigen_vec2
    			right_vec[i][j] = conj( eigen_vec1[i][ dim_i * j + state_id ] );	// the state_id'th col of ConjTranspose[eigen_vec1]

    			left_vec0[i][j] = left_vec[i][j];
    			right_vec0[i][j] = right_vec[i][j];
    		}
    	}
    }
}

double cCLUSTER::calc_Signal_SGL(double * timeseq)
{
	int i, j, k;
	int interval_num, count, count1;
	double sign, e_j, res = 0.0;
	CLD res_amplitude;


	//non-symmetric pulse
	interval_num = pulse_num + 1 ;
	count = 0;	count1 = 0;
	for (i = 0; i < interval_num; i ++)
	{
		count1 >= interval_num ? sign = 1.: sign = -1.;
		for(k = 0; k < blk_num; k ++)
		{
			if(need_calc[k] >= 0)
			{
				for(j = 0; j < blk_dim[k]; j ++)
				{
					i_phase_fact[count] = CLD(0., sign * timeseq[2 * i] * eigen_val1[k][j]);
					count ++;
				}
			}
			else
			{
				for(j = 0; j < blk_dim[k]; j ++)
				{
					i_phase_fact[count] = CLD(0., 0.);
					count ++;
				}
			}
		}
		count1 ++;

		count1 >= interval_num ? sign = 1.: sign = -1.;
		for(k = 0; k < blk_num; k ++)
		{
			if(need_calc[k] >= 0)
			{
				for(j = 0; j < blk_dim[k]; j ++)
				{
					i_phase_fact[count] = CLD(0., sign * timeseq[2 * i + 1] * eigen_val2[k][j]);
					count ++;
				}
			}
			else
			{
				for(j = 0; j < blk_dim[k]; j ++)
				{
					i_phase_fact[count] = CLD(0., 0.);
					count ++;
				}
			}
		}
		count1 ++;
	}
	vzExp(count,i_phase_fact,exp_eigen_val);

	int dim_i;
	for(i = 0; i < blk_num; i++)
    {
    	if(need_calc[i] >= 0)
    	{
    		dim_i = blk_dim[i];

    		for(j = 0; j < dim_i; j++)
    		{
    			left_vec[i][j] = left_vec0[i][j];
    			right_vec[i][j] = right_vec0[i][j];
    		}
    	}
    }

	for(i = 0; i < pulse_num; i ++)
	{
		row_transform_SGL(exp_eigen_val + dim * 2 * i);
		W_Mat_transform_SGL(1);
		row_transform_SGL(exp_eigen_val + dim * (2 * i + 1));
		W_Mat_transform_SGL(2);

	}

	row_transform_SGL(exp_eigen_val + dim * 2 * pulse_num);
	W_Mat_transform_SGL(1);
	row_transform_SGL(exp_eigen_val + dim * (2 * pulse_num + 1));


	res_amplitude = COMPLEX_ZERO;
	for(i = 0; i < blk_num; i++)
	{
		if(need_calc[i] >= 0)
		{
			for(j = 0; j < blk_dim[i]; j ++)
				res_amplitude += left_vec0[i][j] * right_vec[i][j];
		}
	}

	res = abs(res_amplitude);
	return res;
}

void cCLUSTER::row_transform_SGL(CLD * exp_vals)
{

	int i, j, dim_i;
	for(i = 0; i < blk_num; i ++)
	{
		if(need_calc[i] >= 0)
		{
			dim_i = blk_dim[i];

			for(j = 0; j < dim_i; j++)
			{
				right_vec[i][j] *= exp_vals[j];
			}
		}
		exp_vals += blk_dim[i];
	}

}

void cCLUSTER::W_Mat_transform_SGL(int isTrans)
{
	int i, j, dim_i;
	CLD * res;

    char * transa;
	isTrans == 1 ? transa = "NoTranspose" : transa = "ConjugateTranspose";

	for(i = 0; i < blk_num; i++)
	{
		if(need_calc[i] >= 0)
		{
			dim_i = blk_dim[i];
			res = new CLD [dim_i];

			MKL_mat_vec_mult(dim_i, transa, COMPLEX_ONE, W_matrix[i], right_vec[i], COMPLEX_ZERO, res);

			for(j = 0; j < dim_i; j ++)
				right_vec[i][j] = res[j];
		}
	}

	delete [] res;

}

void cCLUSTER::set_localfield_SGL()
{

	int i, j, idx_i, idx_j;
	double z;

	for (i = 0; i < nspin; i++)
	{
		idx_i = spin_index[i];

		z	= spinlist[ idx_i ].Zeeman();

		locX  [i] = spinlist[ idx_i ].mf().x;
		locY  [i] = spinlist[ idx_i ].mf().y;
		locZ  [i] = spinlist[ idx_i ].mf().z + z;
		for (j = 0; j < nspin; j++)
		{
			idx_j = spin_index[j];

			locX  [i] -= zx1[i * nspin + j] * spinlist[ idx_j ].state() * 0.5;
			locY  [i] -= yz1[i * nspin + j] * spinlist[ idx_j ].state() * 0.5;
			locZ  [i] -= zz1[i * nspin + j] * spinlist[ idx_j ].state() * 0.5;
        }
	}

}

//////////////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////////////
void coutHami(double * hami, int dim, int idx)
{
	int i;
	char str[100];

    std::stringstream ss;
    ss << idx;
    ss >> str;

    char dbg_file[200];
    strcpy(dbg_file, WORKSHOP);
    strcat(dbg_file, DEBUG_PATH);
    strcat(dbg_file, "hami");
    strcat(dbg_file, str);
    strcat(dbg_file, ".dat");

    ofstream dbg(dbg_file);
    if(!dbg) assert(0);

    dbg << setprecision(10) <<  setw(15) << setiosflags(ios::scientific);
    for (i = 0; i < dim*dim; i ++)
    {
        dbg << hami[2*i] <<"\t" << hami[2*i+1] << endl;
    }
    dbg << endl;
}

void cout_C_Hami(CLD ** C_HAMI, int blk, int * dim, int idx, int sign)
{
	int i, j, k;
	char str1[100],str2[100];

    std::stringstream ss1, ss2;
    ss1 << idx;		ss2 << sign;
    ss1 >> str1;	ss2 >> str2;

    char dbg_file[200];
    strcpy(dbg_file, WORKSHOP);
    strcat(dbg_file, DEBUG_PATH);
    strcat(dbg_file, "hami");
    strcat(dbg_file, str1);
    strcat(dbg_file, "_");
    strcat(dbg_file, str2);
    strcat(dbg_file, ".dat");
    ofstream dbg(dbg_file);

    if(!dbg) assert(0);
    dbg << setprecision(10) <<  setw(15) << setiosflags(ios::scientific);
    dbg <<"===================================="<< endl;

    for(i = 0; i < blk; i++)
    {
        for (j = 0; j < dim[i]; j ++)
        {
            for(k = 0; k < dim[i]; k++)
                dbg << C_HAMI[i][j*dim[i] + k] << "\t";
            dbg << endl;
        }
        dbg <<"===================================="<< endl;
    }

}

double Lorentzian(double omega, double omega0, double gamma)
{
	double res, dw;
	dw = omega - omega0;
	res = dw * dw + gamma * gamma;
	res = gamma / res;
	return res;
}
//void cout_complex_mat(gsl_matrix_complex * mat, int n1, int n2)
//{
//	for (int i = 0; i < n1; i++)
//	{
//		for (int j = 0; j < n2; j++)
//		{
//			cout << gsl_matrix_complex_get(mat, i, j) << "\t";
//		}
//		cout << endl;
//	}
//}
